pcode)順で各区分を因子化してあります。
オリジナルのデータがどのようになっているかskimrパッケージを用いてサマライズしてみます。元がJSON形式なので、読み込んだ直後は殆どの変量(フィーチャー)が文字型になっているのと欠損値が多いことが分かります。
df %>%
skimr::skim()
| Name | Piped data |
| Number of rows | 88458 |
| Number of columns | 23 |
| _______________________ | |
| Column type frequency: | |
| character | 19 |
| logical | 3 |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| patientId | 0 | 1.00 | 1 | 8 | 0 | 87039 | 0 |
| dateAnnounced | 0 | 1.00 | 10 | 10 | 0 | 253 | 0 |
| gender | 9795 | 0.89 | 1 | 1 | 0 | 2 | 0 |
| detectedPrefecture | 0 | 1.00 | 3 | 15 | 0 | 49 | 0 |
| patientStatus | 84677 | 0.04 | 8 | 23 | 0 | 8 | 0 |
| notes | 46068 | 0.48 | 1 | 270 | 0 | 39760 | 1 |
| mhlwPatientNumber | 88009 | 0.01 | 1 | 11 | 0 | 434 | 0 |
| prefecturePatientNumber | 8324 | 0.91 | 5 | 20 | 0 | 80125 | 0 |
| prefectureSourceURL | 57259 | 0.35 | 5 | 224 | 0 | 3424 | 0 |
| residence | 17671 | 0.80 | 1 | 38 | 0 | 1407 | 0 |
| sourceURL | 586 | 0.99 | 1 | 239 | 0 | 7101 | 0 |
| relatedPatients | 79229 | 0.10 | 2 | 259 | 0 | 5645 | 0 |
| knownCluster | 86146 | 0.03 | 3 | 88 | 0 | 222 | 0 |
| detectedCityTown | 65777 | 0.26 | 2 | 22 | 0 | 660 | 0 |
| cityPrefectureNumber | 65842 | 0.26 | 1 | 34 | 0 | 22607 | 2 |
| citySourceURL | 76727 | 0.13 | 9 | 317 | 0 | 3595 | 0 |
| deceasedDate | 86836 | 0.02 | 10 | 10 | 0 | 203 | 0 |
| deceasedReportedDate | 87358 | 0.01 | 10 | 10 | 0 | 175 | 0 |
| deathSourceURL | 87452 | 0.01 | 14 | 123 | 0 | 618 | 0 |
Variable type: logical
| skim_variable | n_missing | complete_rate | mean | count |
|---|---|---|---|---|
| confirmedPatient | 0 | 1 | 0.98 | TRU: 87038, FAL: 1420 |
| charterFlightPassenger | 88444 | 0 | 1.00 | TRU: 14 |
| cruisePassengerDisembarked | 88447 | 0 | 1.00 | TRU: 11 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| ageBracket | 0 | 1 | 33.54 | 23.17 | -1 | 20 | 30 | 50 | 100 | ▃▇▅▂▁ |
各変量(フィーチャー)を適切な形式に変換し、地域区分でも分析できるように都道府県データと結合します。
x <- df %>%
dplyr::mutate(dateAnnounced = lubridate::as_date(dateAnnounced),
ageBracket = forcats::as_factor(ageBracket),
gender = forcats::as_factor(gender),
patientStatus = forcats::as_factor(patientStatus),
residence = forcats::as_factor(residence),
cluster = dplyr::if_else(!is.na(knownCluster), TRUE, FALSE)) %>%
dplyr::select(dateAnnounced, ageBracket, gender, detectedPrefecture,
patientStatus, knownCluster, cluster) %>%
dplyr::left_join(prefs, by = c("detectedPrefecture" = "pref"))
x
x %>%
skimr::skim()
| Name | Piped data |
| Number of rows | 88458 |
| Number of columns | 13 |
| _______________________ | |
| Column type frequency: | |
| character | 2 |
| Date | 1 |
| factor | 9 |
| logical | 1 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| detectedPrefecture | 0 | 1.00 | 3 | 15 | 0 | 49 | 0 |
| knownCluster | 86146 | 0.03 | 3 | 88 | 0 | 222 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dateAnnounced | 0 | 1 | 2020-01-15 | 2020-10-07 | 2020-08-06 | 253 |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| ageBracket | 0 | 1.00 | FALSE | 13 | 20: 21929, 30: 13552, 40: 11018, -1: 9876 |
| gender | 9795 | 0.89 | FALSE | 2 | M: 44322, F: 34341 |
| patientStatus | 84677 | 0.04 | FALSE | 8 | Dec: 1612, Hos: 1266, Hom: 316, Dis: 283 |
| pcode | 1029 | 0.99 | FALSE | 47 | 13: 27332, 27: 10978, 14: 7351, 23: 5563 |
| 都道府県 | 1029 | 0.99 | FALSE | 47 | 東京都: 27332, 大阪府: 10978, 神奈川: 7351, 愛知県: 5563 |
| 八地方区分 | 1029 | 0.99 | FALSE | 8 | 関東地: 45816, 近畿地: 17587, 九州地: 10027, 中部地: 9012 |
| 広域圏 | 6124 | 0.93 | FALSE | 8 | 首都圏: 46020, 近畿圏: 17045, 中部圏: 7638, 九州圏: 7292 |
| 通俗的区分 | 1029 | 0.99 | FALSE | 11 | 関東: 45816, 関西: 17045, 東海: 7313, 九州: 7292 |
| fct_pref | 1029 | 0.99 | FALSE | 47 | Tok: 27332, Osa: 10978, Kan: 7351, Aic: 5563 |
Variable type: logical
| skim_variable | n_missing | complete_rate | mean | count |
|---|---|---|---|---|
| cluster | 0 | 1 | 0.03 | FAL: 86146, TRU: 2312 |
文字型を因子型に変換するだけでも、変量(フィーチャー)内の水準の比率が大まかにつかめるようになります。
例えば、年代別の陽性判定者数は20代が最も多く、続いて、30代、40代と働き盛りの世代に多いことが分かります。また、都道府県別では、東京、大阪、神奈川、愛知と成っていますが、地方区分別では、以外にも関東、大阪の次に九州がきていることが分かります。
sec_scale <- 50
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced) %>%
dplyr::summarise(n = n()) %>%
dplyr::mutate(n = tidyr::replace_na(n, 0)) %>%
dplyr::mutate(cum = cumsum(n),
ma = zoo::rollmeanr(n, k = 7L, na.pad = TRUE)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced)) +
ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity",
fill = "dark green", alpha = 0.5) +
ggplot2::geom_line(ggplot2::aes(y = cum / sec_scale),
colour = "dark green") +
ggplot2::geom_line(ggplot2::aes(y = ma), colour = "dark red") +
ggplot2::scale_y_continuous(
name = "陽性確定数(日別)・7日間移動平均(濃赤折線)",
sec.axis = ggplot2::sec_axis(~ . * sec_scale,
name = "累積陽性確定数(折線)")
) +
ggplot2::labs(title = paste0("Covid19, Japan(@", Sys.time(), ")"),
subtitle = "", x = "日付")
## Warning: Removed 6 row(s) containing missing values (geom_path).
# plotly::ggplotly()
都道府県別では区分が多すぎるので八地方区分別の積上げ棒グラフを描いてみます。サマリで見たように九州地方が以外にも多いことが読み取れます。
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced, `八地方区分`) %>%
dplyr::summarise(n = n()) %>%
dplyr::mutate(n = tidyr::replace_na(n, 0)) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) +
ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
ggplot2::scale_fill_brewer(palette = "Set1")
折線グラフで表示してみます。
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced, `八地方区分`) %>%
dplyr::summarise(n = n()) %>%
dplyr::mutate(n = tidyr::replace_na(n, 0)) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) +
ggplot2::geom_line(ggplot2::aes(y = n)) +
ggplot2::scale_colour_brewer(palette = "Set1")
さらに分かりやすくするために関東、中部、近畿、九州の四地区を除く他の地区をその他としてまとめてみます。
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
`その他地方` = c("北海道地方",
"東北地方",
"中国地方",
"四国地方"))) %>%
dplyr::group_by(dateAnnounced, `八地方区分`) %>%
dplyr::summarise(n = n()) %>%
dplyr::mutate(n = tidyr::replace_na(n, 0)) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, fill = `八地方区分`)) +
ggplot2::geom_bar(ggplot2::aes(y = n), stat = "identity", alpha = 1.0) +
ggplot2::scale_fill_brewer(palette = "Set1")
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::mutate(`八地方区分` = forcats::fct_collapse(`八地方区分`,
`その他地方` = c("北海道地方",
"東北地方",
"中国地方",
"四国地方"))) %>%
dplyr::group_by(dateAnnounced, `八地方区分`) %>%
dplyr::summarise(n = n()) %>%
dplyr::mutate(n = tidyr::replace_na(n, 0)) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, colour = `八地方区分`)) +
ggplot2::geom_line(ggplot2::aes(y = n)) +
ggplot2::scale_colour_brewer(palette = "Set1") +
ggplot2::facet_wrap(~ `八地方区分`, ncol = 2)
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced, `八地方区分`) %>%
dplyr::summarise(n = n()) %>%
dplyr::ungroup() %>%
tidyr::pivot_wider(names_from = `八地方区分`, values_from = n) %>%
dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>%
tidyr::pivot_longer(cols = -dateAnnounced,
names_to = "八地方区分", values_to = "n") %>%
dplyr::mutate(`八地方区分` = forcats::fct_relevel(`八地方区分`,
"北海道地方", "東北地方",
"関東地方", "中部地方",
"近畿地方", "中国地方",
"四国地方", "九州地方")) %>%
dplyr::group_by(`八地方区分`) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
fill = `八地方区分`)) +
ggplot2::geom_area(stat = "identity", alpha = 1,
position = ggplot2::position_stack(reverse = FALSE)) +
ggplot2::scale_fill_brewer(palette = "Set3") +
ggplot2::labs(title = paste0("Covid19(@", Sys.time(), ")"),
y = "累計人数")
九州が8月頃から急増しているのは、県別に見ると福岡と沖縄での急増が原因と分かります。
x %>%
dplyr::filter(`八地方区分` == "九州地方") %>%
dplyr::group_by(dateAnnounced, `都道府県`) %>%
dplyr::summarise(n = n()) %>%
dplyr::ungroup() %>%
tidyr::pivot_wider(names_from = `都道府県`, values_from = n) %>%
dplyr::mutate_if(is.integer, list(~tidyr::replace_na(., 0L))) %>%
tidyr::pivot_longer(cols = -dateAnnounced,
names_to = "都道府県", values_to = "n") %>%
dplyr::mutate(`都道府県` = forcats::fct_relevel(`都道府県`,
"福岡県", "佐賀県", "長崎県",
"熊本県", "大分県", "宮崎県",
"鹿児島県", "沖縄県")) %>%
dplyr::group_by(`都道府県`) %>%
dplyr::mutate(cum = cumsum(n)) %>%
ggplot2::ggplot(ggplot2::aes(x = dateAnnounced, y = cum,
fill = `都道府県`)) +
ggplot2::geom_area(stat = "identity", alpha = 1,
position = ggplot2::position_stack(reverse = FALSE)) +
ggplot2::scale_fill_brewer(palette = "Set3") +
ggplot2::labs(title = paste0("Covid19, 九州地方(@", Sys.time(), ")"),
y = "累計人数")
x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(`八地方区分`, cluster) %>%
dplyr::summarise(n = n()) %>%
tidyr::pivot_wider(names_from = cluster, values_from = n) %>%
dplyr::mutate(ratio = (`TRUE` / (`TRUE` + `FALSE`) * 100) %>% round(1)) %>%
dplyr::rename(`地方` = `八地方区分`,
`非クラスタ感染者` = `FALSE`, `クラスタ感染者` = `TRUE`,
`クラスタ比率[%]` = ratio)
日毎の陽性判定者数を時系列データ形式に変換します。判定者数がゼロの日は報告されていないので、最初の陽性判定者が出た日から日日次のカレンダーデータを作成して結合してます。なお、時系列データ形式の周期については日次データなので7日間を周期として設定しておきます(考え方によっては約1ヶ月ごとで4週間=28日を周期にするのもありかと思います)。
ts_x <- x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced) %>%
dplyr::summarise(n = n()) %>%
dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>%
dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>%
dplyr::select(-date) %>%
ts(frequency = 1)
tsw_x <- x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced) %>%
dplyr::summarise(n = n()) %>%
dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>%
dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>%
dplyr::select(-date) %>%
ts(frequency = 7)
tsm_x <- x %>%
dplyr::filter(!is.na(fct_pref)) %>%
dplyr::group_by(dateAnnounced) %>%
dplyr::summarise(n = n()) %>%
dplyr::left_join(date, ., by = c("date" = "dateAnnounced")) %>%
dplyr::mutate(n = tidyr::replace_na(n, 0L)) %>%
dplyr::select(-date) %>%
ts(frequency = 28)
時系列データに変換したものをプロットすると可視化の項でプロットした棒グラフと同じ形のグラフになることが分かります。
ts_x %>%
plot()
tsw_x %>%
plot()
tsm_x %>%
plot()
上記からトレンド(長期的傾向)を除いたグラフ。
ts_x %>%
base::diff() %>%
plot()
tsw_x %>%
base::diff() %>%
plot()
tsm_x %>%
base::diff() %>%
plot()
tsw_x %>%
stats::decompose() %>%
plot()
tsm_x %>%
stats::decompose() %>%
plot()
スペクトル密度を推定するためのスペクトル密度の推定
ts_x %>%
spec.pgram()
tsw_x %>%
spec.pgram()
tsm_x %>%
spec.pgram()
自己回帰によるスペクトラム解析
spectrum(ts_x, method="ar")
spectrum(tsw_x, method="ar")
spectrum(tsm_x, method="ar")
ARIMA(自己回帰和分移動平均)モデルによる予測を行ってみます。予測に必要なパラメータはステップワイズにより自動的に最適なものが選択されます。
ts_x %>%
forecast::auto.arima() %>%
forecast::forecast() %>%
plot()
tsw_x %>%
forecast::auto.arima() %>%
forecast::forecast() %>%
plot()
tsm_x %>%
forecast::auto.arima() %>%
forecast::forecast() %>%
plot()